## ============================================================================
## 03_si_s1_measurement.R
## ----------------------------------------------------------------------------
## Purpose:   Assess measurement properties (internal reliability, construct
##            validity via CFA, and inter-item correlations) for all
##            multi-item scales in Study 1 and Study 2.
##            Corresponds to SI Sections S1.3.1 and S1.4.1.
##
## Inputs:    combined_data.rds (via 00_setup.R)
## Outputs:   Figures/fig_s4.pdf,
##            Figures/fig_s5.pdf,
##            Tables/table_s6.tex,
##            Tables/table_s7.tex,
##            Tables/table_s8.tex,
##            Tables/table_s9.tex,
##            Tables/table_s10.tex
## ============================================================================

source("00_setup.R")

#### ============== Study 1: Audience Effects ============== ####
## ---- Table S6: Internal reliability by sample -------
alpha_dat <- list()

## Collective Guilt
alpha_guilt_au <-
  combined_dat %>%
  filter(sample == "AU", Z_land_info == 0) %>%
  select(y_guilt_q1_n, y_guilt_q2_n, y_guilt_q3_n)

pcor <- psych::polychoric(alpha_guilt_au)$rho
omega <- MBESS::ci.reliability(alpha_guilt_au, type = "categorical")

alpha_dat[[1]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_guilt_au)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_guilt_au))$total[1]),
  omega = omega$est,
  outcome = "Collective Guilt",
  sample = "Australia"
)

alpha_guilt_us <-
  combined_dat %>%
  filter(sample == "US", Z_land_info == 0) %>%
  select(y_guilt_q1_n, y_guilt_q2_n, y_guilt_q3_n)

pcor <- psych::polychoric(alpha_guilt_us)$rho
omega <- MBESS::ci.reliability(alpha_guilt_us, type = "categorical")

alpha_dat[[2]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_guilt_us)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_guilt_us))$total[1]),
  omega = omega$est,
  outcome = "Collective Guilt",
  sample = "United States"
)

## Colorblind Racism
alpha_cbr_au <-
  combined_dat %>%
  filter(sample == "AU", Z_land_info == 0) %>%
  select(y_cbr_q1_n, y_cbr_q2_n, y_cbr_q3_n)

pcor <- psych::polychoric(alpha_cbr_au)$rho
omega <- MBESS::ci.reliability(alpha_cbr_au, type = "categorical")

alpha_dat[[3]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_cbr_au)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_cbr_au))$total[1]),
  omega = omega$est,
  outcome = "Colorblind Racism",
  sample = "Australia"
)

alpha_cbr_us <-
  combined_dat %>%
  filter(sample == "US", Z_land_info == 0) %>%
  select(y_cbr_q1_n, y_cbr_q2_n, y_cbr_q3_n)

pcor <- psych::polychoric(alpha_cbr_us)$rho
omega <- MBESS::ci.reliability(alpha_cbr_us, type = "categorical")

alpha_dat[[4]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_cbr_us)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_cbr_us))$total[1]),
  omega = omega$est,
  outcome = "Colorblind Racism",
  sample = "United States"
)

## Symbolic Reparations
alpha_symbolic_au <-
  combined_dat %>%
  filter(sample == "AU", Z_land_info == 0) %>%
  select(y_symbolic_q1_n, y_symbolic_q2_n, y_symbolic_q3_n,
         y_symbolic_q4_n)

pcor <- psych::polychoric(alpha_symbolic_au)$rho
omega <- MBESS::ci.reliability(alpha_symbolic_au, type = "categorical")

alpha_dat[[5]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_symbolic_au)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_symbolic_au))$total[1]),
  omega = omega$est,
  outcome = "Symbolic Reparations",
  sample = "Australia"
)

alpha_symbolic_us <-
  combined_dat %>%
  filter(sample == "US", Z_land_info == 0) %>%
  select(y_symbolic_q1_n, y_symbolic_q2_n, y_symbolic_q3_n,
         y_symbolic_q4_n)

pcor <- psych::polychoric(alpha_symbolic_us)$rho
omega <- MBESS::ci.reliability(alpha_symbolic_us, type = "categorical")

alpha_dat[[6]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_symbolic_us)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_symbolic_us))$total[1]),
  omega = omega$est,
  outcome = "Symbolic Reparations",
  sample = "United States"
)

## Substantive Reparations
alpha_substantive_au <-
  combined_dat %>%
  filter(sample == "AU", Z_land_info == 0) %>%
  select(y_substantive_q1_n, y_substantive_q2_n, y_substantive_q3_n,
         y_substantive_q4_n)

pcor <- psych::polychoric(alpha_substantive_au)$rho
omega <- MBESS::ci.reliability(alpha_substantive_au, type = "categorical")

alpha_dat[[7]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_substantive_au)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_substantive_au))$total[1]),
  omega = omega$est,
  outcome = "Substantive Reparations",
  sample = "Australia"
)

alpha_substantive_us <-
  combined_dat %>%
  filter(sample == "US", Z_land_info == 0) %>%
  select(y_substantive_q1_n, y_substantive_q2_n, y_substantive_q3_n,
         y_substantive_q4_n)

pcor <- psych::polychoric(alpha_substantive_us)$rho
omega <- MBESS::ci.reliability(alpha_substantive_us, type = "categorical")

alpha_dat[[8]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_substantive_us)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_substantive_us))$total[1]),
  omega = omega$est,
  outcome = "Substantive Reparations",
  sample = "United States"
)

## Export results
bind_rows(alpha_dat) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  select(sample, outcome, alpha, palpha, omega) %>%
  pivot_wider(names_from = c(sample), values_from = c(alpha, palpha, omega)) %>%
  select(outcome, ends_with("_Australia"), everything()) %>%
  xtable() %>%
  print(include.rownames = FALSE, include.colnames = FALSE, only.contents = TRUE,
        type = "latex", sanitize.text.function = identity, comment = FALSE,
        hline.after = NULL, file = paste0(tab_dir, "table_s6.tex"))


## ---- Figure S4: Correlations among individual items by sample ------
# Define items and labels
item_vars <- c(
  "y_guilt_q1_n", "y_guilt_q2_n", "y_guilt_q3_n",
  "y_cbr_q1_n", "y_cbr_q2_n", "y_cbr_q3_n",
  "y_symbolic_q1_n", "y_symbolic_q2_n", "y_symbolic_q3_n", "y_symbolic_q4_n",
  "y_substantive_q1_n", "y_substantive_q2_n", "y_substantive_q3_n", "y_substantive_q4_n"
)

item_labels <- c(
  paste0("Guilt-", 1:3),
  paste0("CBR-", 1:3),
  paste0("SymRep-", 1:4),
  paste0("SubRep-", 1:4)
)

# Function to compute polychoric matrix and convert to tidy format
get_polychoric_long <- function(data, sample_name) {
  rho_mat <- psych::polychoric(data)$rho
  as.data.frame(rho_mat) %>%
    rownames_to_column("Var1") %>%
    pivot_longer(-Var1, names_to = "Var2", values_to = "Correlation") %>%
    mutate(
      Var1 = factor(Var1, levels = item_vars, labels = item_labels),
      Var2 = factor(Var2, levels = item_vars, labels = item_labels),
      row_num = as.numeric(Var1),
      col_num = as.numeric(Var2)
    ) %>%
    filter(row_num > col_num) %>%
    mutate(sample = sample_name)
}

# Create AU and US datasets subsetting to untreated respondents
au_items <- combined_dat %>%
  filter(sample == "AU", Z_land_info == 0) %>%
  select(all_of(item_vars))

us_items <- combined_dat %>%
  filter(sample == "US", Z_land_info == 0) %>%
  select(all_of(item_vars))

# Get tidy correlation matrices
au_long <- get_polychoric_long(au_items, "Australia")
us_long <- get_polychoric_long(us_items, "United States")

# Combine
rho_long_all <- bind_rows(au_long, us_long)

# Plot
p <-
  ggplot(rho_long_all, aes(x = Var1, y = Var2, fill = Correlation)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", Correlation)), size = 3, color = "black") +
  scale_fill_gradient2(
    low = "blue", high = "red", mid = "white", midpoint = 0,
    limit = c(-1, 1), name = "Corr"
  ) +
  facet_col(~ sample) +
  theme_tufte(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank(),
    strip.text = element_text(size = 14, hjust = 0.5, face = "bold"),
    plot.margin = unit(c(
      t = 0.05,
      r = 0.5,
      b = 0.05,
      l = 0.05
    ), "lines"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank(),
  ) +
  labs(x = NULL, y = NULL)

ggsave(
  filename = paste0(fig_dir, "fig_s4.pdf"),
  p,
  width = 6.5,
  height = 6
)

## ---- In-text: Sample sizes for Figure S4 caption ----
nrow(au_items) # Australia untreated n
nrow(us_items) # United States untreated n

## ---- Confirmatory factor analysis by sample -------

us_data <-
  combined_dat %>%
  filter(sample == "US", Z_land_info == 0)

au_data <-
  combined_dat %>%
  filter(sample == "AU", Z_land_info == 0)

## Collective guilt:
guilt_model <- 'collective_guilt =~ y_guilt_q1_n + y_guilt_q2_n +  y_guilt_q3_n'

guilt_fit_us <- cfa(guilt_model, data = us_data, std.lv = TRUE, ordered = TRUE)

guilt_fit_au <- cfa(guilt_model, data = au_data, std.lv = TRUE, ordered = TRUE)

## NB: 3-items -> saturated model so fit stats are not meaningful

# Extract standardized solution
guilt_std_us <- standardizedSolution(guilt_fit_us)
guilt_std_au <- standardizedSolution(guilt_fit_au)

## Colorblind racism:
cbr_model <- 'colorblind_racism =~ y_cbr_q1_n + y_cbr_q2_n + y_cbr_q3_n'

cbr_fit_us <- cfa(cbr_model, data = us_data, std.lv = TRUE, ordered = TRUE)

cbr_fit_au <- cfa(cbr_model, data = au_data, std.lv = TRUE, ordered = TRUE)

# Extract standardized solution
cbr_std_us <- standardizedSolution(cbr_fit_us)
cbr_std_au <- standardizedSolution(cbr_fit_au)

## Symbolic reparations:
symbolic_model <- 'symbolic =~ y_symbolic_q1_n + y_symbolic_q2_n + y_symbolic_q3_n + y_symbolic_q4_n'

symbolic_fit_us <- cfa(symbolic_model, data = us_data, std.lv = TRUE,
                       ordered = TRUE)

symbolic_fit_au <- cfa(symbolic_model, data = au_data, std.lv = TRUE,
                       ordered = TRUE)

# Extract standardized solution
symbolic_std_us <- standardizedSolution(symbolic_fit_us)
symbolic_std_au <- standardizedSolution(symbolic_fit_au)

## Substantive reparations:
substantive_model <- 'substantive =~ y_substantive_q1_n + y_substantive_q2_n + y_substantive_q3_n + y_substantive_q4_n'

substantive_fit_us <- cfa(substantive_model, data = us_data, std.lv = TRUE,
                          ordered = TRUE, estimator = "WLSMV")

substantive_fit_au <- cfa(substantive_model, data = au_data, std.lv = TRUE,
                          ordered = TRUE, estimator = "WLSMV")

# Extract standardized solution
substantive_std_us <- standardizedSolution(substantive_fit_us)
substantive_std_au <- standardizedSolution(substantive_fit_au)

cfa_dat <-
  bind_rows(
    guilt_std_au %>% mutate(sample = "AU", y = "Collective Guilt"),
    guilt_std_us  %>% mutate(sample = "US", y = "Collective Guilt"),
    cbr_std_au %>% mutate(sample = "AU", y = "Colorblind Racism"),
    cbr_std_us  %>% mutate(sample = "US", y = "Colorblind Racism"),
    symbolic_std_au %>% mutate(sample = "AU", y = "Symbolic Reparations"),
    symbolic_std_us  %>% mutate(sample = "US", y = "Symbolic Reparations"),
    substantive_std_au %>% mutate(sample = "AU", y = "Substantive Reparations"),
    substantive_std_us  %>% mutate(sample = "US", y = "Substantive Reparations")
  ) %>%
  mutate(
    statistic = case_when(
      op == "=~" ~ "Loading",
      op == "~~" &
        lhs == rhs & grepl("y_", lhs) ~ "Residual Variance"
    ),
    sample = factor(
      sample,
      levels = c("AU", "US"),
      labels = c("Australia", "United States")
    ),
    y = factor(
      y,
      levels = c(
        "Collective Guilt",
        "Colorblind Racism",
        "Symbolic Reparations",
        "Substantive Reparations"
      )
    ),
    item = rhs,
    item_short = str_replace(item, "^y_(guilt|cbr|symbolic|substantive)_q(\\d+)_n$", function(x) {
      m <- str_match(x, "^y_(guilt|cbr|symbolic|substantive)_q(\\d+)_n$")
      prefix <- case_match(
        m[, 2],
        "symbolic" ~ "sym",
        "substantive" ~ "sub",
        .default = m[, 2]
      )
      paste0(prefix, "_", m[, 3])
    })
  ) %>%
  filter(statistic %in% c("Loading", "Residual Variance")) %>%
  select(sample,
         y,
         item,
         item_short,
         statistic,
         est.std,
         se,
         ci.lower,
         ci.upper) 

## NB: Residual Variance = 1 - Loading^2

## ---- In-text: CFA fit statistics ----
fitMeasures(substantive_fit_us, c("cfi", "tli", "srmr"))
fitMeasures(substantive_fit_au, c("cfi", "tli", "srmr"))

fitMeasures(symbolic_fit_us, c("cfi", "tli", "srmr"))
fitMeasures(symbolic_fit_au, c("cfi", "tli", "srmr"))
## ---- Tables S7-S8: Factor loadings by sample -----

## Collective Guilt and Colorblind Racism (3-item scales)
cfa_dat %>%
  filter(statistic == "Loading",
         !str_detect(y, "Reparations")) %>%
  mutate(
    prefix = str_extract(item_short, "^[^_]+"),
    item_num = str_extract(item_short, "\\d+$"),
    item_label = paste0("Item ", item_num),
  ) %>%
  select(sample, y, item_label, est.std) %>%
  pivot_wider(names_from = c(sample, y), values_from = est.std) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  select(item_label, starts_with("Australia_"), everything()) %>%
  xtable() %>%
  print(include.rownames = FALSE, include.colnames = FALSE, only.contents = TRUE,
        type = "latex", sanitize.text.function = identity, comment = FALSE,
        hline.after = NULL, file = paste0(tab_dir, "table_s7.tex"))

## Substantive/Symbolic Reparations (4-item scales)
cfa_dat %>%
  filter(statistic == "Loading",
         str_detect(y, "Reparations")) %>%
  mutate(
    prefix = str_extract(item_short, "^[^_]+"),
    item_num = str_extract(item_short, "\\d+$"),
    item_label = paste0("Item ", item_num),
  ) %>%
  select(sample, y, item_label, est.std) %>%
  pivot_wider(names_from = c(sample, y), values_from = est.std) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  select(item_label, starts_with("Australia_"), everything()) %>%
  xtable() %>%
  print(include.rownames = FALSE, include.colnames = FALSE, only.contents = TRUE,
        type = "latex", sanitize.text.function = identity, comment = FALSE,
        hline.after = NULL, file = paste0(tab_dir, "table_s8.tex"))



#### ============== Study 2: Practitioner Effects ============== ####
## Setup: construct stacked vignette data
mine_dat <-
  combined_dat %>%
  select(
    sample,
    admin_response_id,
    y_trust_m_n,
    y_reputation_m_n,
    y_effective_m_n,
    y_sincerity_m_n,
    Z_vignette_a,
    Z_vignette_order,
    Z_vig_first,
    Z_vig_second,
    Z_vig_first_type,
    Z_vig_second_type
  ) %>%
  group_by(sample) %>%
  mutate(
    y_mine_aidx = c(
      y_trust_m_n + y_reputation_m_n + y_effective_m_n +
        y_sincerity_m_n
    ),
    y_mine_aidx_s = glass_delta(Y = y_mine_aidx, Z = Z_vignette_a,
                                reference = "Control"),
    Z_vig_type = "Mining"
  ) %>%
  rename(y_aidx = y_mine_aidx,
         y_aidx_s = y_mine_aidx_s,
         Z = Z_vignette_a,
         Z_vig_order = Z_vignette_order,
         y_trust_n = y_trust_m_n,
         y_reputation_n = y_reputation_m_n,
         y_effective_n = y_effective_m_n,
         y_sincerity_n = y_sincerity_m_n)

ent_dat <-
  combined_dat %>%
  select(
    sample,
    admin_response_id,
    y_trust_e_n,
    y_reputation_e_n,
    y_effective_e_n,
    y_sincerity_e_n,
    Z_vignette_b,
    Z_vignette_order,
    Z_vig_first,
    Z_vig_second,
    Z_vig_first_type,
    Z_vig_second_type
  ) %>%
  group_by(sample) %>%
  mutate(
    y_ent_aidx = c(
      y_trust_e_n + y_reputation_e_n + y_effective_e_n +
        y_sincerity_e_n
    ),
    y_ent_aidx_s = glass_delta(Y = y_ent_aidx, Z = Z_vignette_b,
                               reference = "Control"),
    Z_vig_type = "Entertainment"
  ) %>%
  rename(y_aidx = y_ent_aidx,
         y_aidx_s = y_ent_aidx_s,
         Z = Z_vignette_b,
         Z_vig_order = Z_vignette_order,
         y_trust_n = y_trust_e_n,
         y_reputation_n = y_reputation_e_n,
         y_effective_n = y_effective_e_n,
         y_sincerity_n = y_sincerity_e_n
  )

stack_dat <-
  bind_rows(mine_dat, ent_dat) %>%
  group_by(sample) %>%
  ## Standardize for pooled estimation (i.e., intercept is control mean)
  mutate(y_aidx_s2 = glass_delta(Y = y_aidx, Z = Z,
                                 reference = "Control")) %>%
  ungroup()

## ---- Table S9: Internal reliability by sample -------
alpha_dat <- list()

alpha_csr_au <-
  stack_dat %>%
  ungroup() %>%
  filter(Z == "Control", sample == "AU") %>%
  select(y_trust_n, y_reputation_n, y_effective_n, y_sincerity_n)

pcor <- psych::polychoric(alpha_csr_au)$rho
omega <- MBESS::ci.reliability(alpha_csr_au, type = "categorical")

alpha_dat[[1]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_csr_au)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_csr_au))$total[1]),
  omega = omega$est,
  outcome = "Corporate Social Responsibility",
  sample = "Australia"
)

alpha_csr_us <-
  stack_dat %>%
  ungroup() %>%
  filter(Z == "Control", sample == "US") %>%
  select(y_trust_n, y_reputation_n, y_effective_n, y_sincerity_n)

pcor <- psych::polychoric(alpha_csr_us)$rho
omega <- MBESS::ci.reliability(alpha_csr_us, type = "categorical")

alpha_dat[[2]] <- data.frame(
  alpha = as.numeric(psych::alpha(alpha_csr_us)$total[1]),
  palpha = as.numeric(psych::alpha(pcor, n.obs = nrow(alpha_csr_us))$total[1]),
  omega = omega$est,
  outcome = "Corporate Social Responsibility",
  sample = "United States"
)

## Export results
bind_rows(alpha_dat) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  select(sample, alpha, palpha, omega) %>%
  xtable() %>%
  print(include.rownames = FALSE, include.colnames = FALSE, only.contents = TRUE,
        type = "latex", sanitize.text.function = identity, comment = FALSE,
        hline.after = NULL, file = paste0(tab_dir, "table_s9.tex"))

## ---- Figure S5: Correlations among individual items by sample ------
# Define items and labels
item_vars <- c(
 "y_trust_n", "y_reputation_n", "y_effective_n", "y_sincerity_n"
)

item_labels <- paste0("CSR-", 1:4)

# Function to compute polychoric matrix and convert to tidy format
get_polychoric_long <- function(data, sample_name) {
  rho_mat <- psych::polychoric(data)$rho
  as.data.frame(rho_mat) %>%
    rownames_to_column("Var1") %>%
    pivot_longer(-Var1, names_to = "Var2", values_to = "Correlation") %>%
    mutate(
      Var1 = factor(Var1, levels = item_vars, labels = item_labels),
      Var2 = factor(Var2, levels = item_vars, labels = item_labels),
      row_num = as.numeric(Var1),
      col_num = as.numeric(Var2)
    ) %>%
    filter(row_num > col_num) %>%
    mutate(sample = sample_name)
}

# Create AU and US datasets subsetting to control group obs.
au_items <-
  stack_dat %>%
  ungroup() %>%
  filter(Z == "Control", sample == "AU") %>%
  select(all_of(item_vars))

us_items <-
  stack_dat %>%
  ungroup() %>%
  filter(Z == "Control", sample == "US") %>%
  select(all_of(item_vars))

# Get tidy correlation matrices
au_long <- get_polychoric_long(au_items, "Australia")
us_long <- get_polychoric_long(us_items, "United States")

# Combine
rho_long_all <- bind_rows(au_long, us_long)

# Plot
p <-
  ggplot(rho_long_all, aes(x = Var1, y = Var2, fill = Correlation)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", Correlation)), size = 3, color = "black") +
  scale_fill_gradient2(
    low = "blue", high = "red", mid = "white", midpoint = 0,
    limit = c(-1, 1), name = "Corr"
  ) +
  facet_row(~ sample) +
  theme_tufte(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank(),
    strip.text = element_text(size = 14, hjust = 0.5, face = "bold"),
    plot.margin = unit(c(
      t = 0.05,
      r = 0.5,
      b = 0.05,
      l = 0.05
    ), "lines"),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank(),
  ) +
  labs(x = NULL, y = NULL)

ggsave(
  filename = paste0(fig_dir, "fig_s5.pdf"),
  p,
  width = 5,
  height = 3
)

## ---- Confirmatory factor analysis by sample -------

us_data <-
  stack_dat %>%
  filter(sample == "US", Z == "Control")

au_data <-
  stack_dat %>%
  filter(sample == "AU", Z == "Control")

## Collective guilt:
csr_model <- 'y_csr =~ y_trust_n + y_reputation_n + y_effective_n + y_sincerity_n'

csr_fit_us <- cfa(csr_model, data = us_data, std.lv = TRUE, ordered = TRUE)

csr_fit_au <- cfa(csr_model, data = au_data, std.lv = TRUE, ordered = TRUE)

# Extract standardized solution
csr_std_us <- standardizedSolution(csr_fit_us)
csr_std_au <- standardizedSolution(csr_fit_au)

cfa_dat <-
  bind_rows(
    csr_std_au %>% mutate(sample = "AU", y = "Corporate Social Responsibility"),
    csr_std_us  %>% mutate(sample = "US", y = "Corporate Social Responsibility")
  ) %>%
  mutate(
    statistic = case_when(
      op == "=~" ~ "Loading",
      op == "~~" &
        lhs == rhs & grepl("y_", lhs) ~ "Residual Variance"
    ),
    sample = factor(
      sample,
      levels = c("AU", "US"),
      labels = c("Australia", "United States")
    ),
    item = rhs
  ) %>%
  filter(statistic %in% c("Loading", "Residual Variance")) %>%
  select(sample,
         y,
         item,
         statistic,
         est.std,
         se,
         ci.lower,
         ci.upper) %>%
  mutate(across(where(is.numeric), \(x) round(x, 3)))

## NB: Residual Variance = 1 - Loading^2

## ---- In-text: CFA fit statistics ----
fitMeasures(csr_fit_us, c("cfi", "tli", "srmr"))
fitMeasures(csr_fit_au, c("cfi", "tli", "srmr"))

## ---- Table S10: Factor loadings by sample -----

cfa_dat %>%
  filter(statistic == "Loading") %>%
  mutate(item_label = paste0("Item ", match(item, unique(item)))) %>%
  select(sample, y, item_label, est.std) %>%
  pivot_wider(names_from = c(sample, y), values_from = est.std) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  select(item_label, starts_with("Australia_"), everything()) %>%
  xtable() %>%
  print(include.rownames = FALSE, include.colnames = FALSE, only.contents = TRUE,
        type = "latex", sanitize.text.function = identity, comment = FALSE,
        hline.after = NULL, file = paste0(tab_dir, "table_s10.tex"))




## ---- Session Info -----------------------------------------------------------
sessionInfo()